{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "<center>Автор материала: Владимир Яшин (@vdyashin), студент-магистр ВШЭ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center> Метод прыжков для выбора числа кластеров"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Известно, что, помимо обучения с учителем, в машинном обучении есть обучение без учителя. Задачи кластеризации являются задачами без учителя. Есть ряд методов кластеризации, которые требуют указать число кластеров (см. Spectral clustering, Ward hierarchical clustering, K-Means, Agglomerative clustering и другие методы). А что если мы изначально не знаем количество кластеров?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Примером области, где важно узнать число кластеров, может быть криминология. Криминалистам может понадобиться узнать количество людей, которые причастны к написанию как-то письма, будь-то электроенного или нет.  Урбанистам выбор числа кластеров может помочь понять, сколько разных районов в некоторой агломерации. Возможно, даже предсказать, в каких округах может возникнуть гетто. Такую задачу часто решают маркетологи, когда им нужно узнать, на сколько кластеров поделить своих клиентов, чтобы проводить таргетированные рекламные акции."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "На одном из занятий научно-исследовательского семинара в университете мне было предложено реализовать один из алгоритмов, которые позволяют определить число кластеров в датасете. Алгоритм был предложен в статье:\n",
    "\n",
    "    Sugar, C. A., & James, G. M. (2003). Finding the number of clusters in a dataset: An information-theoretic \n",
    "    approach. Journal of the American Statistical Association, 98(463), 750–763. JOUR.\n",
    "\n",
    "Метод назвается методом \"прыжков\". Если коротко описывать алгоритм работы, то мы для выбранного диапазона $K$ в алгоритме $K$-средних, считаем внитриклассовую дисперсию, а потом смотрим ее прирост, подсчитанный с некоторым гипер-параметорм, и выбираем там, где такой прирост (\"прыжок\") был наибольший. Если ничего не понятно, то прочитайте этот небольшой мануал и вы, скорее всего, разберетесь."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Статья, которую вы сейчас читаете состоит из практической части, в которой метод будет использован скорее как \"черный ящик\". Прочитав только практическую часть, вы получите некоторое поверхностное представление о методе. Это ваша \"программа минимум\". Если захочется большего, то теоретическая часть, которая следует за практической, даст вам понимание того, как устроен алгоритм. Для того, чтобы понять _почему_ он работает, прочтайте исходную статью, ее можно найти в моем репозитории https://github.com/vdyashin/JumpMethod/blob/master/article/article.pdf."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Практика"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Давайте начнем с простых примеров, где мы знаем количество кластеров и распределения, из которых данные будем сэмплировать. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сейчас мы создадим четыре кластера на плоскости с известными нам средними и их матрицей ковариации."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Начнем с параметров. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# means\n",
    "Mu1, Mu2, Mu3, Mu4 = [-4, -4], [-4, 4], [4, -4], [4, 4]\n",
    "\n",
    "# covariance\n",
    "r = 0\n",
    "\n",
    "# covariance matrix\n",
    "Sigma = [[1, r], \n",
    "         [r, 1]]\n",
    "\n",
    "# inverse the cov matrix\n",
    "Sigma_inv = np.linalg.inv(Sigma)\n",
    "\n",
    "# sample size\n",
    "n = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь насэмплируем наблюдения. Тут нам поможет `numpy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample\n",
    "c1 = np.random.multivariate_normal(mean = Mu1, cov = Sigma, size = n)\n",
    "c2 = np.random.multivariate_normal(mean = Mu2, cov = Sigma, size = n)\n",
    "c3 = np.random.multivariate_normal(mean = Mu3, cov = Sigma, size = n)\n",
    "c4 = np.random.multivariate_normal(mean = Mu4, cov = Sigma, size = n)\n",
    "\n",
    "# concatenate them into one dataset\n",
    "C = np.row_stack([c1, c2, c3, c4])\n",
    "print(C.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Давайте посмотрим на получившийся датасет. Используем `matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(C[:, 0], C[:, 1], 'o', alpha=0.2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В этом случае очень просто определить количество кластеров, даже не зная, что их было четыре. Однако, давайте проверим работу нашего алгоритма."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сейчас я не хочу, чтобы вы вдавались в подробности того, как алгоритм работает. Если заинтересует, то почитайте потом раздел [теории](#theory), а пока просто пролистните следующую ячейку, ну или просто запустите, если открыли в ноутбуке. <a id='class'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import numpy as np\n",
    "# from sklearn.cluster import KMeans\n",
    "\n",
    "class JumpsMethod(object):\n",
    "    \"\"\"\n",
    "    This class is initialized with dataset (data) and has two methods that help\n",
    "    to calculate distortions and jumps -- the concepts proposed in \n",
    "        Sugar, C. A., & James, G. M. (2003). Finding the number of clusters in \n",
    "        a dataset: An information-theoretic approach. Journal of the American \n",
    "        Statistical Association, 98(463), 750–763. JOUR.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, data):\n",
    "        self.data = data\n",
    "        # dimension of 'data'; data.shape[0] would be size of 'data'\n",
    "        self.p = data.shape[1]\n",
    "        # vector of variances (1 by p)\n",
    "        \"\"\" 'using squared error rather than Mahalanobis distance' (SJ, p. 12)\n",
    "        sigmas = np.var(data, axis=0)\n",
    "        ## by following the authors we assume 0 covariance between p variables (SJ, p. 12)\n",
    "        # start with zero-matrix (p by p)\n",
    "        self.Sigma = np.zeros((self.p, self.p), dtype=np.float32)\n",
    "        # fill the main diagonal with variances for\n",
    "        np.fill_diagonal(self.Sigma, val=sigmas)\n",
    "        # calculate the inversed matrix\n",
    "        self.Sigma_inv = np.linalg.inv(self.Sigma)\"\"\"\n",
    "    \n",
    "    \n",
    "    def distortions(self, cluster_range=(1, 10), random_state=0):\n",
    "        \"\"\" returns a vector of calculated distortions for each cluster number.\n",
    "            If the number of clusters is 0, distortion is 0 (SJ, p. 2) \n",
    "            'cluster_range' -- range of numbers of clusters for KMeans;\n",
    "            'data' -- n by p array \"\"\"\n",
    "        # dummy vector for Distortions\n",
    "        self.distortions = np.repeat(0, cluster_range[1] + 1).astype(np.float32)\n",
    "\n",
    "        # for each k in cluster range implement\n",
    "        for k in range(cluster_range[0], cluster_range[1] + 1):\n",
    "            # initialize and fit the clusterer giving k in the loop\n",
    "            KM = KMeans(n_clusters=k, random_state=random_state)\n",
    "            KM.fit(self.data)\n",
    "            # calculate centers of suggested k clusters\n",
    "            centers = KM.cluster_centers_\n",
    "            # since we need to calculate the mean of mins create dummy vec\n",
    "            for_mean = np.repeat(0, len(self.data)).astype(np.float32)\n",
    "\n",
    "            # for each observation (i) in data implement\n",
    "            for i in range(len(self.data)):\n",
    "                # dummy for vec of distances between i-th obs and k-center\n",
    "                dists = np.repeat(0, k).astype(np.float32)\n",
    "\n",
    "                # for each cluster in KMean clusters implement\n",
    "                for cluster in range(k):\n",
    "                    # calculate the within cluster dispersion\n",
    "                    tmp = np.transpose(self.data[i] - centers[cluster])\n",
    "                    \"\"\" 'using squared error rather than Mahalanobis distance' (SJ, p. 12)\n",
    "                    dists[cluster] = tmp.dot(self.Sigma_inv).dot(tmp)\"\"\"\n",
    "                    dists[cluster] = tmp.dot(tmp)\n",
    "\n",
    "                # take the lowest distance to a class\n",
    "                for_mean[i] = min(dists)\n",
    "\n",
    "            # take the mean for mins for each observation\n",
    "            self.distortions[k] = np.mean(for_mean) / self.p\n",
    "\n",
    "        return self.distortions\n",
    "    \n",
    "    \n",
    "    def jumps(self, Y=None):\n",
    "        \"\"\" returns a vector of jumps for each cluster \"\"\"\n",
    "        # if Y is not specified use the one that suggested by the authors (SJ, p. 2) \n",
    "        if Y is None:\n",
    "            self.Y = self.p / 2\n",
    "        \n",
    "        else:\n",
    "            self.Y = Y\n",
    "        \n",
    "        # the first (by convention it is 0) and the second elements\n",
    "        self.jumps = [0] + [self.distortions[1] ** (-self.Y) - 0]\n",
    "        self.jumps += [self.distortions[k] ** (-self.Y) \\\n",
    "                       - self.distortions[k-1] ** (-self.Y) \\\n",
    "                       for k in range(2, len(self.distortions))]\n",
    "        \n",
    "        # calculate recommended number of clusters\n",
    "        self.recommended_cluster_number = np.argmax(np.array(self.jumps))\n",
    "        \n",
    "        return self.jumps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Начнем с инициализации класса. Класс принимает на вход данные."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "JM = JumpsMethod(data=C)\n",
    "print('Количество признаков:', JM.p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Далее, прежде чем переходить к подсчету \"прыжков\", которые позволят нам показать, сколько кластеров оптимально выбрать, нужно подсчитать отклонения (distortions). В данном случае мы переберем с одного до десяти кластеров. Подсчитав в классе отклонения, можно посмотреть на их график в зависимости от числа кластеров."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "JM.distortions(cluster_range=(1, 10), random_state=13);\n",
    "\n",
    "plt.plot(JM.distortions, '-o')\n",
    "plt.ylabel('Distortion')\n",
    "plt.xlabel('Number of Centers')\n",
    "plt.xlim(1, 10)\n",
    "plt.ylim(0, 20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Наконец, подсчитаем значения \"прыжков\" и посмотрим на график."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "JM.jumps()\n",
    "\n",
    "plt.plot(JM.jumps, '-o')\n",
    "plt.ylabel('Jump')\n",
    "plt.xlabel('Number of Centers')\n",
    "plt.xlim(1, 10)\n",
    "plt.ylim(0, 1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видим, что такой пример не вызывает затруднений у алгоритма. Давайте теперь попробуем более сложный пример. Будем использовать те же параметры, кроме средних значений."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# means\n",
    "Mu1, Mu2, Mu3, Mu4 = [-2, -2], [-2, 2], [2, -2], [2, 2]\n",
    "\n",
    "# covariance\n",
    "r = 0\n",
    "\n",
    "# covariance matrix (Gamma)\n",
    "Sigma = [[1, r], \n",
    "         [r, 1]]\n",
    "# inversed cov matrix\n",
    "Sigma_inv = np.linalg.inv(Sigma)\n",
    "\n",
    "# sample size\n",
    "n = 1000\n",
    "\n",
    "c1 = np.random.multivariate_normal(mean = Mu1, cov = Sigma, size = n)\n",
    "c2 = np.random.multivariate_normal(mean = Mu2, cov = Sigma, size = n)\n",
    "c3 = np.random.multivariate_normal(mean = Mu3, cov = Sigma, size = n)\n",
    "c4 = np.random.multivariate_normal(mean = Mu4, cov = Sigma, size = n)\n",
    "\n",
    "C = np.row_stack([c1, c2, c3, c4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(C[:, 0], C[:, 1], 'o', alpha=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "JM = JumpsMethod(data=C)\n",
    "\n",
    "JM.distortions(cluster_range=(1, 10), random_state=13);\n",
    "\n",
    "plt.plot(JM.distortions, '-o')\n",
    "plt.ylabel('Distortion')\n",
    "plt.xlabel('Number of Centers')\n",
    "plt.xlim(1, 10)\n",
    "plt.ylim(0, 20)\n",
    "plt.show();\n",
    "\n",
    "JM.jumps()\n",
    "\n",
    "plt.plot(JM.jumps, '-o')\n",
    "plt.ylabel('Jump')\n",
    "plt.xlabel('Number of Centers')\n",
    "plt.xlim(1, 10)\n",
    "plt.ylim(0, 1);\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Уже менее уверенно, но все же точно. Теперь пример еще сложнее."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# means\n",
    "Mu1, Mu2, Mu3, Mu4 = [-1, -1], [-1, 1], [1, -1], [1, 1]\n",
    "\n",
    "# covariance\n",
    "r = 0\n",
    "\n",
    "# covariance matrix (Gamma)\n",
    "Sigma = [[1, r], \n",
    "         [r, 1]]\n",
    "# inversed cov matrix\n",
    "Sigma_inv = np.linalg.inv(Sigma)\n",
    "\n",
    "# sample size\n",
    "n = 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c1 = np.random.multivariate_normal(mean = Mu1, cov = Sigma, size = n)\n",
    "c2 = np.random.multivariate_normal(mean = Mu2, cov = Sigma, size = n)\n",
    "c3 = np.random.multivariate_normal(mean = Mu3, cov = Sigma, size = n)\n",
    "c4 = np.random.multivariate_normal(mean = Mu4, cov = Sigma, size = n)\n",
    "\n",
    "C = np.row_stack([c1, c2, c3, c4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(C[:, 0], C[:, 1], 'o', alpha=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Границы кластера вовсе не различимы. Посмотрим, что покажут результаты алгоритма."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "JM = JumpsMethod(data=C)\n",
    "\n",
    "JM.distortions(cluster_range=(1, 10), random_state=13);\n",
    "\n",
    "plt.plot(JM.distortions, '-o')\n",
    "plt.ylabel('Distortion')\n",
    "plt.xlabel('Number of Centers')\n",
    "plt.xlim(1, 10)\n",
    "plt.ylim(0, 20)\n",
    "plt.show();\n",
    "\n",
    "JM.jumps()\n",
    "\n",
    "plt.plot(JM.jumps, '-o')\n",
    "plt.ylabel('Jump')\n",
    "plt.xlabel('Number of Centers')\n",
    "plt.xlim(1, 10)\n",
    "plt.ylim(0, 1);\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из последнего графика видно, что алгоритм совсем не справляется и показывает, что лучше всего принять количество кластеров за единицу. Чисто визуально, там действительно один кластер."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Теория <a id='theory'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для того, чтобы перейти к обсуждению алгоритма, необходимо сперва дать определение наиболее важному и, как мне кажется, наиболее трудному для восприятия концепту отклонений (distortions), которые, по сути, являются внутриклассовой дисперсией. \n",
    "\n",
    "Положим $ D $ наши данные состоят из $ p $ признаков, каждый из которых набран из распредения с матрицей ковариации $ \\Sigma $. Тогда $ x_i \\in D $ многомерная случайная величина размерности $ p $. Так же определим $ c_1, c_2, \\dots, c_K $ как центры кластеров для соответствующей случайной величины, причем, $ c_{x_i} $ ближайший центр кластера для $ x_i $. Тогда наименьшее достижимое отклонение для $ K $, полученных при применении метода K-средних\n",
    "\n",
    "$$\n",
    "d_K = \\frac{1}{p} \\min_{c_1, c_2, \\dots, c_K} \\mathbb{E} \\Big[\\big(x_i - c_{x_i}\\big)^T \\Sigma^{-1} \\big(x_i - c_{x_i}\\big)\\Big]\n",
    "$$\n",
    "\n",
    "что является средним расстоянием Махаланобиса между $ x_i $ и $ c_{x_i} $ по всем $ p $ признакам. Очевидно, что нам редко известна матрица ковариаций $ \\Sigma $, поэтому авторы метода предлагают использовать квадратичную ошибку (p. 12).\n",
    "\n",
    "Обратите внимание, что на самом деле нужно подсчитать среднее минимумов, а не минимум средних. Возможно авторы просто ошиблись. Ну или я, но по-другому не работает."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь перейдем непосредственно к алгоритму подсчета \"прыжков\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Необходимо подсчитать $ \\hat{d}_K $ для каждого из рассматриваемого количества кластеров, то есть нужно запустить метод $K$-средних для любого на ваш выбор диапазона $K$. На выходе получится вектор из отклонений длины, равной количеству перебираемых значений кластеров.\n",
    "2. Выбрать гипер-параметр $ Y $, степень трансформации (transformation power). Авторы советуют начать подбор с $ Y = \\frac{p}{2}$.\n",
    "3. Подсчитать \"прыжки\" по формуле: $ J_K = \\hat{d}^{-Y}_{K} - \\hat{d}^{-Y}_{K-1} $. Проще говоря, это просто прирост отклонений. Будьте внимательны, $ d_0 := 0$.\n",
    "4. Принять за $K^*$ такое количество кластеров, при котором $J$ максимальное. Технически: $ K^* = \\arg \\max_K J_K $. Обратите внимание, что метод может выбрать $ K^* = 1 $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Моя попытка написать класс, который исполняет этот алгоритм, показана [выше](#class).\n",
    "\n",
    "На этом все! Спасибо, если прочитал не только Практику, но и [Теорию](#theory). \n",
    "\n",
    "В этой статье я показал, как можно выбирать количество кластеров для алгоритма K-средних. Более детальное объяснение метода смотрите в исходной статье. Если захочется что-то законтрибьютить или поправить зайдите сюда: https://github.com/vdyashin/JumpMethod."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
